ог 


УДК 004.93'1;004.932 


Д.А. Гончаров, А.М. Недзьведь, С.В. Абламейко 
Объединенный институт проблем информатики НАН РЬБ, г. Минск, Беларусь 
гопсвагоуда(@ета!.сот, педтуеда@пеултап.Баз-пе.Бу 


Анализ методов построения объемных 
* 
дистанционных карт 


В статье рассмотрены и сопоставлены друг с другом по нескольким аспектам (ТОЧНОСТЬ, скорость, простота 
реализации, возможность распараллеливания) методы построения трехмерных дистанционных карт, в 
результате чего можно с определенной долей уверенности дать ответ на вопрос, в каком случае лучше 
использовать тот или иной метод. 


Введение 


Дистанционная карта (дистанционное преобразование) — это такое представле- 
ние данных, в котором для каждой точки объекта известно расстояние до ближайшей 
точки фона. Впервые в литературе по обработке изображений дистанционные преоб- 
разования появились в 1966 году в работе Розенфельда и Пфалца [1], в которой они 
применили дистанционное преобразование Чамфера для построения скелета объекта. 
С тех пор был не только улучшен тот алгоритм, но и были придуманы альтернативные 
решения. Были разработаны алгоритмы для различных представлений данных: шести- 
гранных сеток, точечных решеток, продолговатых сеток. Дистанционные преобразо- 
вания применяются во многих областях обработки изображений (морфологическая 
обработка [2], скелетонизация и построение срединных осей [3]), визуализации (уско- 
рение метода бегущего луча, эндоскопия), моделировании (обнаружение столкновений, 
морфинг, анимация, конструктивная блочная геометрия). В связи с увеличивающимся 
применением трехмерных объектов в различных областях науки, возникла необходи- 
мость в быстрых и точных алгоритмах для построения объемных дистанционных карт. 

Целью данной работы является систематизация известных алгоритмов постро- 
ения дистанционных карт, выявление наиболее быстрых и точных из них, а также 
оценка области дальнейшего применения каждого из них. 


Классификация алгоритмов 


Алгоритмы для построения карт расстояний могут быть классифицированы по 
нескольким критериям. 

Классификация по способу вычисления расстояния: 
— Дистанционные преобразования Чамфера, где новое значение расстояния для вок- 
селя вычисляется из расстояний его соседей и соответствующих весов маски. 
— Векторные дистанционные преобразования, в которых для каждого обработан- 
ного вокселя хранится вектор ближайших точек поверхности, а для необработанного 
вокселя этот вектор строится, используя вектора соседей из шаблонной окрестности. 


"Работы выполнялись в рамках проекта ТЗТС В-1486. 
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— Дистанционные преобразования на основе решения уравнения эйконала, где рас- 
стояние находится из конечных разностей первого или второго порядков расстояний 
соседних вокселей. 
— Евклидово дистанционное преобразование в квадрате, в котором вместо рассто- 
яния до ближайшей точки фона хранится квадрат этого расстояния, что позволяет 
воспользоваться некоторыми преимуществами. 

По схеме обхода классификация выглядит следующим образом: 
— Последовательные, в которых обход начинается из одного угла и идет в направлении 
противоположного угла. Обычно требуют несколько проходов по всему объему данных. 
— Волновые, в которых вычисление расстояний начинается от границы в порядке уве- 
личения расстояния. Если требуется, вычисления можно остановить на определенном 
значении расстояния. Это свойство волновой схемы дает ряд преимуществ, например, 
для вычисления морфологических операторов, построения множеств уровня. 


Дистанционные преобразования Чамфера 


В дистанционных преобразованиях Чамфера (ДПЧ) [4], [5] шаблон (рис. 1) цент- 
рирован относительно каждого вокселя и расстояние вычисляется как минимум среди 
сумм расстояний соседних вокселей и соответствующего шаблонного веса. Для преоб- 
разований этого вида разработаны обе схемы обхода (последовательная и волновая). 


Е Е О ес а 
ГЕ е Е ес се ав а 6 а 
19) фе] ав а а а 0 
ГЕ е Е есь |с 
Е Е Геае Е 
7-2 7-1 7 


Рисунок 1 — Шаблон расстояний для ДП Чамфера для прямого обхода 


В последовательной схеме используются два прохода: прямой и обратный. На 
первом проходе ищутся минимальные расстояния до нижней границы изображения, 
на обратном - до верхней. Далее результаты обоих проходов объединяются и полу- 
чается результат. 


Таблица 1 — Веса, используемые для разных метрик и размеров структурного элемента 
Метрика, размер структурного а Ь с (| е Е 
элемента 
Городская 1 
Шахматная 1 1 
Евклидова, 3 Ч 33 
Квази-Евклидова, 5° те 3: ||:5/5 6 3 


В волновой схеме воксели обрабатываются в порядке увеличения расстояния. 
Чтобы гарантировать правильный порядок обработки, обработанные воксели поме- 
щаются в приоритетную очередь. На каждом шаге достается воксель с головы очереди и 
для всех необработанных соседних вокселей вычисляется расстояние по заданному 
шаблону и они вставляются в хвост очереди. Алгоритм заканчивает работу, когда 
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очередь пуста. Асимптотическая сложность такой схемы порядка О(М№10$ №) ‚ что хуже, 
чем в последовательной схеме (О(М№)). Однако, при специальных ограничениях на ве- 


личину необходимого расстояния, такой подход оправдывает себя и работает значительно 
быстрее, чем схема с просмотром всех вокселей. 


Векторные дистанционные преобразования 


ДП Чамфера обладают невысокой точностью, что особенно проявляется на больших 
расстояниях от поверхности объекта. Эту проблему помогают решить векторные ДП 
(ВДП) [6], в которых для каждого вокселя хранится вектор ближайших точек границы. 

Рассмотрим ВДП [7]. Изначально вокселям объекта приписывается расстояние 0, 
остальным — бесконечность. Всего требуется 8 обходов по всему объему данных (4 пря- 
мых и 4 обратных) (рис. 2). При каждом обходе обновляется вектор ближайших точек 
границы, используя вектора соседей. 

Существуют также другие шаблоны и схемы обхода для разных метрик. Некото- 
рые позволяют делать обходы независимо друг от друга, что дает возможность распа- 
раллелить процесс вычисления карты расстояний. 
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Рисунок 2 — Шаблоны для векторных ДП для прямого (слева) 
и обратного (справа) проходов 


Быстрый маршевый метод и его вариации 


Пусть у нас имеется множество ©, ОО -— его граница. Уравнением эйконала 
называется уравнение в частных производных вида: 
УТО] Е) = ЕО 


) 1 
ТхебдО =0 — 


где У означает градиент, | — какую-то метрику, Е(х) — скорость движения точки (функ- 


ция с положительными значениями. Т(х) называется временем прибытия и имеет 
следующий физический смысл: это минимальное время, за которое граница достигает 
точки х. Если положить Ё(х) =1, то значение Т(х) будет не что иное, как расстояние 
до ближайшей точки границы множества. 

Идея алгоритма [8] заключается в том, что мы расширяем границу с постоян- 
ной единичной скоростью в направлении нормали к ней. По своей сути это метод с 


волновой схемой обхода и его структура почти идентична алгоритму Дейкстры нахож- 
дения кратчайшего пути. В классическом алгоритме быстрой маршировки расстояния 
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пересчитываются с помощью формулы конечных разностей первого порядка. Для на- 
хождения вокселя с минимальным значением расстояния используется приоритетная 
очередь. Сложность управления приоритетной очередью порядка О(М№1юз №), где М- 
число вокселей. Как отмечалось выше, в целом этот алгоритм, как и все алгоритмы с 
волновой схемой, работает медленнее, чем последовательные ДП. 


Рисунок 3 — Пример работы быстрого маршевого метода. Изначально 
просматриваются все воксели, смежные с граничными (черные клетки). Затем 
выбирается клетка с наименьшим значением расстояния и помечается 
как просмотренная 


Также для пересчета расстояний можно использовать формулу конечных раз- 
ностей второго порядка, что увеличивает точность вычислений. 

Яцив и др. [9] заменили кучу массивом связанных списков. Операции добавле- 
ния и извлечения выполняются за О(1), но из-за того, что списки несортированные, 
возникают некоторые погрешности. В некоторых случаях они приемлемы, учитывая 
экономию времени. Но на сложных объектах ошибки становятся существенными. 


Евклидово дистанционное преобразование в квадрате 


Этот алгоритм был разработан Мейстером, Рординком и Хесселинком [10]. Ос- 
новная идея — использование дистанционных карт меньших размерностей. 
Опишем алгоритм для двухмерного изображения. Трехмерный случай аналоги- 
чен. 
По определению: 
й а и 2 о 
1а(х, у) =шт{(х-П+(-Л) у, ЛЕО, (2) 
где /[(1, /) — пиксель с координатами (2, /). 
Формула (2) может быть переписана следующим образом: 
р д2 у 2 
1а(х, у) = пп (х-й +0), (3) 


где 


8 у) =пу- Л, ЦЬЛ ЕО. (4) 

Алгоритм состоит из двух фаз. Первая фаза использует 2 прохода для вычис- 
ления дистанционного преобразования по столбцам (функция ©). Это несложно 
сделать с помощью динамического программирования за линейное время. 

Во второй фазе идет вычисление строка за строкой для всех х с фиксированной 
у -координатой. Поэтому можем ввести следующее обозначение ©(1) = С(р, у). 

Затем: 

4(ь у) = шт /65, (5) 


где 


Л(хл) = (х-0’ +80. (6) 
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Полезно рассмотреть геометрическую интерпретацию формулы (6). График 
Л(х,?) — это парабола с вершиной ([,2(1)). И (5) мы можем интерпретировать как 


нижнюю границу семейства парабол /(х,1) (рис. 4). 


НХ 


Рисунок 4 — Нижняя граница семейства парабол 


Чтобы найти нижнюю границу, все параболы помещаются в стек. Когда рас- 
сматривается очередная парабола, она выталкивает параболы из стека, которые ее выше. 
Дополнительно хранится информация, на каком участке произошло выталкивание. За- 
тем эта парабола заносится в стек. В итоге легко можно определить нижнюю границу. 

Так как вычисления для столбцов (строк) независимы, алгоритм может быть рас- 
параллелен. 


Сравнение 


Алгоритмы сравнивались на нескольких медицинских изображениях. Первый 
тест выполнялся на объектах простой формы: точка в центре области и лимфоузлы, 
полученные с изображений компьютерной томографии (КТ). Второй — для изобра- 
жений сложной формы: кость с полостью (КТ), щитовидная железа (ультразвуковое 
изображение, УЗИ), опухоль забрюшинного пространства (КТ). Третий — для сложных 
ветвящихся объектов: сгенерированное дерево с множеством сегментов небольшой 
толщины, ангиографическое изображение сосудов (КТ), клеточный конгломерат, полу- 
ченный методами кофокальной микроскопии. 

Максимальная ошибка (МО) вычислялась по формуле: 


МЕ(Г) = тах(`/(х) —1 (х)) УхеЕГ, (7) 


где / — построенная карта расстояний, / — эталонная карта. 


еоп 


еоп 
Средняя ошибка вычислялась следующим образом: 


11-1) 


АЕ = (8) 


где № -— общее число вокселей. 
Тестирование проводилось на двуядерном процессоре Соте 2 Рио 1.66 ГГц. 


мух 
а 


Рисунок 5 — Тестовое изображение: кость с полостью внутри и дерево 
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се = в 
Рисунок 6 — Слои тестовых медицинских изображений: компьютерной томографии, 
ультразвукового исследования 
Результаты работы алгоритмов представлены в табл. 2. 


Таблица 2 — Сравнительная таблица алгоритмов (СО — средняя ошибка, 
МО -— максимальная ошибка, В — время выполнения) 


Иа: 1 тест 2 тест 3 тест 
р Со | мо В со | мо | в | СО [| м В 
ВДИ 0 0 7.10 0.003 | 0.097 9.10 0.02 0.31 9.35 


ДИЧ (городская м-ка) 63.11 | 149.7 5.52 | 17.38 | 93.32 | 7.74 | 11.30 | 76.93 7.76 
ДПЧ (@пахматная м-ка) 15.93 | 49.58 6.91 10.39 | 35.79 | 8.67 | 4.16 | 22.38 8.35 


ДПЧ (евклидова м-ка, 33) | 8.48 | 15.03 | 11.70 | 2.51 | 15.84 | 12.33 | 2.04 | 13.14 | 10.85 
ДПЧ (евклидова м-ка, 53) | 2.15 | 6.31 | 35.56 | 0.444 | 3.69 | 48.28 | 0.59 | 4.50 | 3274 


Упрощенный БММ 0.37 0.55 125.8 0.04 0.27 116.2 | 0.06 0.93 130.23 
Упрощенный БММ, т =3 0.01 0.01 0.22 0.026 0.01 3.32 0.05 0.055 3.48 
ДИЧ 
3 0.21 0.41 0.21 0.27 0.65 3.48 0.29 0.71 3/1 
(евклидова м-ка, 3”), т =3 
ЕДП? 0 0 5.80 0 0 8.17 0 0 9.05 


Некоторые результаты сравнений представлены на рис. 7 и рис. 8. 


И Вы ты альт 
вДП ЧДП(3х3х3) Ч1П(5х5х5) МЕМ ЕДП2 


"точка (и) == лимфоузел = кость “щитовидная жедеза 
= ‘опухоль дерево (и) '' сосуды 


Рисунок 7 — Сравнение алгоритмов — средняя ошибка 


150 


ока 


[осы ае ии 


вдП ЧДП(3х3х3) ЧДП(5х5х5) МБМ ЕДП2 


"точка (и) =" лимфоузел =="кость =“ щитовидная жедеза 
= ‘опухоль дерево (и) '' сосуды 


Рисунок 8 — Сравнение алгоритмов -— скорость работы (с) 


На некоторых тестах самыми быстрыми оказались ДП Чамфера. Они также прос- 
ты в реализации, время их работы предсказуемо. ДПЧ с маской размерами 5х5х5 
оказалось точнее ДПЧ с маской размерами 3х3х3, но время в 3 - 4 раза выше. ДИЧ 
могут применяться в областях, где важнее скорость, чем точность (например, при ус- 
корении метода бегущего луча). 

Векторное ДП также было быстрым на некоторых тестах, хотя его точность на- 
много выше, чем у ДП Чамфера. Тем не менее, этот метод не является абсолютно 
точным -— на большом удалении от границы объекта начинают появляться ошибки. 
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Также из-за необходимости хранения векторов ближайших точек алгоритм исполь- 
зует много памяти. Он может применяться в областях, где важны и скорость, и точ- 
ность (обработка столкновений, морфинг). 

Упрощенный быстрый маршевый метод из-за своей асимптотической сложнос- 
ти значительно проигрывает остальным методам, но все же обладает очень хорошей 
точностью на несложных объектах. Благодаря волновой схеме обхода данных, его можно 
использовать там, где интересуемое расстояние ограничено (например, морфологичес- 
кие операторы) и не требуется полный проход по данным. В этих случаях он быстрее 
алгоритмов с последовательной схемой прохода. 


Выводы 


Самым точным преобразованием для всех классов объектов оказалось ЕДП2. 
Однако из-за своей сложности этот алгоритм оказался не самым быстрым -— на 
объектах первой и цельных объектов второй группы (щитовидная железа и опухоли) 
быстрее оказалось ВДП и на более сложных объектах второй группы (кости и опухо- 
ли сложной формы, сосуды) — ЧДП. ЧДП с маской размерами 3х3х3 было быстрее, 
чем ЧДП с маской 5х5х5, но построенная карта расстояний имеет существенно худ- 
шее качество. Результаты работы МБМ самые медленные, но благодаря его волновой 
схеме обхода этот метод можно использовать для анализа очень сложных объектов, 
таких как сосудистые и нервные сети. Для исследования небольших фрагментов изо- 
бражений время работы алгоритма удовлетворительное и позволяет включать его в 
автоматизированные системы морфометрического анализа. 
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Д.А. Гончаров, А.М. Недзьведь, С.В. Абламейко 

Анал метод1в побудови об’есмних дистанщИних карт 

У стати розглянуп 1 зставлен! один 3 одним за деюлькома аспектами (точесть, швидюсть, простота реализаций, 
можлив!сть розпаралелювання) методи побудови тривимрних дистанщйних карт, в результат! чого можна з 
певною долею впевненост! дати в1дпов1дь на питання, в якому випадку краще використовувати той або 
1нший метод. 


р. Сопсйагоъ, А. Мейсуеа, 5. АМатеуко 

Апа1у$15 Мео4$ Гог Сопзбгис_то Уоштейлс О1$апсе Мар5 

Мефод$ оЁ Ф$апсе тар сопзиасноп Уеге ехатше ап согтпраге4 у сасК офег ш зеуега! азресз (ассигасу, 
зрее4, еазе оЁ парететаноп, ве роз Шу оЁ рагаПеИпо) {есблаиез Юг сопзёласйие @гее-ипепзюпа! Чапсе. 
МИ зоте сопйадепсе уоц сап э1уе ап апз\уег ‘0 а даезНоп, пп УЛусь сазе 1 15 Бейег 10 зе опе ог апофег тле#оа. 
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